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Abstract 

We suggest a modification of the operator exponential method for the numerical solv- 
•^T ' ing the difference linear initial boundary value problems. The scheme is based on the 

^^ , representation of the difference operator for given boundary conditions as the per- 

turbation of the same operator for periodic ones. We analyze the error, stability and 
efficiency of the scheme for a model example of the one-dimensional operator of second 



^ 



o 



C^ ' difference. 

1 Introduction 

ff^ , Numerical solution of the linear difference initial-boundary value problems is an essential 

^^ ' part of modelling of the physical processes and phenomena described by the evolutionary 

r^^ ■ differential equations, such as the Schrodinger equation, diffusion equation, Ginzburg- 

^^ , Landau equation and many other. 

Along with the classical grid methods |]l|, an ever increasing use in treatment of such 
■~V ■ evolutionary problems is currently made of the operator exponential (OE) method [^, 

"^ which is based on the Lee-Trotter-Kato formula [^ for approximate calculation of the 

exponential of the sum of noncommuting matrices. 

The OE method offers a number of advantages relevant to both explicit and implicit 

difference schemes. It does not involve iteration procedures and often proves to be ab- 

* ^^ I 

KA , solutely stable. Its applicability is, however, limited by the impossibility to explicitly 

^j ' calculate the exponential of the difference operators expressed in a general form. In fact, 

effective algorithms of exponential calculation exist only for the difference operators with 
constant coefficients and periodic boundary conditions on "rectangular" subsets of Z". 
These algorithms are based on the fast Fourier transform Q| and allow one to calculate 
the exponential in 0{N log2 A^) operations, where N is the number of points in the domain. 
For other boundary conditions such algorithms are not available. 

In this work a linear difference operator with assigned boundary conditions is consid- 
ered as perturbation of the same operator with periodic boundary conditions, and the 
exponential of such an operator is calculated by the Lee-Trotter-Kato formula. 
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The perturbating operator is essentially an operator in the space of functions on the 
domain's boundary, so the problem of calculating its exponential is essentially simpler 
because the number of boundary points is generally much smaller than the total number 
of points in the domain. This ensures practically the same efficiency of the algorithm 
proposed as that obtained with the Fourier method for solving the periodic boundary 
value problems. 

Analysis of the error, stability and efficiency of the algorithm proposed is generally quite 
complicated. So we only present it for a model example of the operator of second difference 
(one-dimensional difference Laplace operator). This case is probably least "favorable" 
for the OE method due to availability of effective difference schemes such as the sweep 
method |jl]. Nevertheless, the algorithm proposed is competitive with the well-known 
schemes, in particular, as applied to the Schrodinger equation. 

The idea of representing differential operators with various boundary conditions as one 
another's perturbations was put forward by M G Krein [^ and is being actively used in 
modern mathematical physics (see, for example, [0]). Applicability of the Krein method 
to difference operators was considered in [^]. The results of this work were reported at 
the "Conference on differential equations and applications" (Saransk, Russia, 1994). The 
summary of this report was published in ||5| . 



2 Description of the method 

Let A^ be a set of points and C{M.) a set of complex-valued functions on M.. Let A : 
C{M) — > C{M) be a linear operator of the form 

(i/)(x) = ^ a,iy)f{y), f G C(A^), (2.1) 

where 'jx is a finite subset of A4 for each value of x, and ax{-) is a given function on jx- 

Let il be a subset of M. We call point x G fi an inner point of 0, relative to A if 
Jx ^ O, and a boundary point of 17 relative to A if jx does not completely lie in 0,. Denote 
by Oa^ the set of all boundary points of 0, relative to A and let 5^0 = U 7x\^- 



Note that by definition (2J.), to calculate the values of Af at boundary points of O, 



we have to know the values of function / on the set U ^yiO. A linear operator 

L : C{n) — >C{nu bA^) such that (L/)(x) = f{x) for all x G J7 

will be called an extension operator for A. An operator Al : (7(^2) — > (7(17) such that 
(^2,/)(x) = {ALf){x) for all X G r^ will be called an L-expansion of operator A. The 
operator L plays the same role for difference operators as the boundary conditions play 
for differential operators. 

We now consider a difference initial-boundary value problem for operator A: 

df 

-L = Aif, t > 0, 

dt '^■'' - ' (2.2) 

/(0,x)=5(x), 



Operator Exponential Method for Initial-Boundary Value Problems 315 

where f{t,-),g G C{^), L is a given extension operator for A. The solution of problem 
dU) is of the form: 

fit, •) = exj>{tAL)g, (2.3) 

where operator exp{tAL) can be defined as a matrix power series since Q is finite. Given A 
and O, the efficiency of computation of exp(i^i) in ( |2.3| ) may largely depend on the 
extension operator L. Let us clarify the above said with a test example. 
Let A4 = Z, let A be the difference Laplacian j^] defined by the relation 

(A/)(x)=/(x-l)-2/(x) + /(x + l), xGZ. (2.4) 

In this case jx in ( ^ ) is the set {x — l,x,x + 1}, 

{1 if y = X — 1, 
-2 if y = X, 
1 if y = x + l. 

Let J7 = {0, 1,... ,iV-l}. ThendA^ = {0,N-l},andbA^ = {-l,N}. Let the extension 
operator L correspond to the periodic boundary conditions for A: 



{Lf){x) 



r fiN- 


-1) if x = -l 


f{x) 


if xen, 


[ m 


if x = N. 



(2.5) 



The exponential exp(tAj^) in (|2.3D can be expressed by the following formula: 

exp(tAL) = P-^ exp(tA)F, (2.6) 

where A is the diagonal operator of the form: 

(A/)(x) = iy{x)f{x), ij{x) = -4sin2 ^, x e Q, 



F is the operator of the discrete Fourier transform: 
(F/)(x) = 5:exp(-^)/(y). 



yen 

Note that computation of the vector exp(tAj^)/ via formula ( |2.6| ) takes about A^log2 N 
operations if we make use of the known Fast Fourier Transform (FFT) algorithm [Q]. 

Let K be the extension operator for A, corresponding to the boundary conditions of 
the 3rd kind, i.e., 

r a/(0) if X = -1, 

iKf)ix) = { fix) if xen, (2.7) 

[ Pf{N - 1) if x = N, 

where a and /? are, generally, the complex coefficients. (The case a = (3 = —1 corresponds 
to the Dirichlet boundary conditions, and a = /3 = 1 to the Neumann boundary condi- 
tions.) In this case the known algorithms for exact computation of the vector exp(tAi^)/ 
(for example, using expansion in eigenfunctions of A^^) involve ~ N'^ operations. 
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Considering the general case again, the question arises: whether the available effective 
algorithm for the exp{tAL) computation (L is the given extension operator) can be used 
to approximately evaluate exp(t^;^) for another extension operator K? 

Below we describe a version of an OE method which establishes the relation between 
the exponents of different extensions of a difference operator and thus answer the above 
question. 

Let K and L be two different extension operators for the operator A. We further 
assume for simplicity that these operators satisfy the following additional condition: the 
equations Lf = Lg and Kf = Kg are fulfilled for any f,gG C{Q) such that f{x) = g{x) 
for X G Q\dA^- 

Consider operator Gkl = Ak — A^. It follows from definition of extension operators, 
that {GKLf){x) = at all inner points x € fi, and that GklJ = Gkl9 if /(x) = g{x) 
at the inner points x G O. Therefore, Gkl is the direct sum of the zero operator in the 
subspace C(0\9a^) of C(0) and an operator in the subspace C(5a^); we will denote the 
restriction of Gkl on (7(9^17) by the same character Gkl- 

This suggests that when the number of boundary points of Jl is much smaller than the 
total number of points in Q, the problem of computing exp(tG'i^L) becomes much simpler 
than the initial problem of evaluating exp(tj4;^). 

Remark. For the extension operators of the general form a small modification of these 
arguments leads to the same result. In the above example of operator A the number of 
boundary points is equal to two, and the computation reduces to finding the exponential 
of a 2 X 2 matrix. 

Since Ak = Al + Gkl, the following relations hold: 

exp(tii^) = exp(tiL) eMtGKL) + O {t^) = Si{t) + O [t^] (2.8) 

exp(tii^) = exp ( -IGkl j exp(tiL) exp i ^tGKL j + O (t^) 

= S2it) + O [t^] . (2.9) 

They are similar to the conventionally used OE schemes of the 1st and 2nd order approxi- 
mation [^. 

Owing to the above mentioned properties of operator Gkl, these formulas allow one to 
roughly calculate the exponential of Ak with almost same efficiency as that of exp(iA^). 
The natural domain of application of ( |2.8D , ( ^.g| ) is the one when A is a difference operator 
with constant coefficients in C(Z^), s > 1 (i.e., the functions ax{y) in expression ( p.l| ) only 

s 

depend on the difference x — y), and ft = f^ {0, 1, . . . , A'j- — 1} is a parallelepiped in Z^. In 

this case there is a specific extension operator L which is defined by the relation {Lf){x) = 
f{x mod A^), where (x mod N)j = Xj mod Nj for j = 0,1, . . . ,s, corresponding to the 
periodic boundary conditions for A. The exact value of operator exp{tAK) is calculated 
using a multidimensional discrete Fourier transform, the calculation procedure involves 
about M log2 M operations, where M = Ni- . . .■ Ns- The number of points of the set d^^ 
can be estimated as 



\dAn\<c{A)MJ2w. 



A,- 
1 ^ 



Operator Exponential Method for Initial-Boundary Value Problems 



317 



where the constant C{A) depends on #(supp a(x)). Hence, if C{A) <^ uiinNj, then 
the size of matrix Gkl is much smaller than that of Ak', this allows to effectively use 
formulas ( |2.8D , (|2.9|). Such a situation occurs in approximations of differential operators 
with difference ones, and the constant C{A) in this case depends on the order of the 
operator approximated and, generally, on a method of approximation. 



3 Error and stability of the method 



Let us now study error and stability of numerical algorithms based on formulas (2.^ 
and (127 



Let f{t, •) = exj){tAK)g be the exact solution of problem (|2.2[) with operator Ak and 
hj{t, •) = Sj{t)g, j = 1,2, g £ C{Q). As an error estimate of one step of the OE algorithms 
we consider the norms of differences of the functions f{t, •) and hj{t, •): 

ej{t,g) = \\f{t,-)-h,{t,-)\\, i = l,2, geC{n). 

By expanding / — hj in the Taylor series at t = we find 



si{t,g) 
£2{t,g) 



Ak,Gkl 



g 



\ + o {t^) 



Ak, 



Ak,Gkl 



Gkl, 



Ak,Gkl 



^ + o('') 



12 



(3.1) 
(3.2) 



where [A, B] = AB - BA. 

Thus, the schemes ( |2.8D and ( |2.9D have the first and second order of approximation, 
respectively. 

For comparison, consider similar estimates for the classical Euler and Krank-Nickolson 
(KN) methods O. The corresponding approximations are of the form (cf. ( p.8| ) and ( |2.9D ) 

exp(iix) = E + tAK + 0{t^) ^ Se{t) + 0{t'^) 



for the Euler scheme and 

exp(iix) = iE + hAK){E - ^tAx)-^ + 0{t^) = Skn{t) + ©(t^) 



(3.3) 



(3.4) 



for the Krank-Nicolson one. It is easy to obtain the estimates for errors of these approxi- 
mation, namely 



t' 



S,{t,g) = \\AU\j + 0{t') 



(3.5) 
(3.6) 



for Euler and KN schemes, respectively. Clearly, unlike the values 6j, the estimates £j, 
where j = 1,2, are determined by the norms of commutators of Ak with Gkl, rather 
than by the powers of Ak- This accounts for the differences in the features of the OE 
algorithms and classical schemes. 
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Stability analysis of the OE methods (12.81) and (|2.9|) requires evaluation of the norms 
of the relevant step operators Si{t) and 5*2 (t). We are going to show that these methods 
are stable for rather small t > 0, if operator Ax satisfies the condition 



Re{AKg, g) < for any g G C{n). 
Indeed, the functions 

Sj{t,g) = \\Sjgf, g£C{n), for j = 1,2 
are analytic in t in a vicinity of zero, and Sj{0,g) = \\g\\'^- Besides, 



(3.7) 



dt 



2Re(ix<?,5)<0 



t=o 



due to the assumption (3.7). Therefore, for a sufficiently small positive values of t the 
inequality 



, „2 dSi 



+ om<\\gr-. 



t=0 



Sj{t,g) 



is valid, i.e., the schemes (|2.8D, ( |2.9D are stable. Note that the condition ( |3.7| ) means that 
the spectrum spec Ak lies in the left half-plane [^ ; this guarantees stability of the initial 
problem (p.2|). 

There exist two classes of operators Ak for which the OE method shows absolute 
stability. 



1. Schemes ( ^.8]) , ( ^.9]) are absolutely stable, if Ak and Al are Hermitian operators 
and spec Al and spec Gkl both lie in the left half-plane. This immediately follows from 
the simplest estimates: 



\Si 

\S2 



< ||exp(tGKL)||||exp(ML)|| < 1, 

2 



< 



exp ( -tCxL 



\exp{tAL)\\ < 1. 



This class of operators includes, in particular, the difference Laplace operator with the 
Dirichlet boundary conditions. 

2. The OE method is also absolutely stable, if operators Ak and Al are both skew- 
symmetric, i.e., Ak = —^*k s-'^d ^l = ~^l- This condition is equivalent to the case 



when specj4i^ and specAi both lie on the imaginary axis. Then operators (2^), (2^) are 
unitary (just as the operator exp{tAK))', hence, their norms are equal to 1. An example 
of such case is the Schrodinger operator. 

The OE algorithm may, however, lack absolute stability even when both spec Ak and 
spec Al are in the left half-plane. This loss of stability is associated with the positive 
eigenvalues available for the "boundary" operator Gkl- The example is a difference 
Laplacian with the Neumann boundary conditions. 
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The above estimates can be illustrated by the Z)-expansion of Laplace operator A^) with 
the extension operator D corresponding to the Dirichlet boundary conditions (a = /? = —1 

are absolutely stable. 



m 



In this case, as mentioned earlier, schemes (2 



The error estimates for the time step in the OE methods ( p.lD and (|3.2| ) depend on the 
initial vector g. As typical vectors g we consider the eigenfunctions of the operator A/j: 




".w-v/^™'""^" 



N 



'4 



for j,k = 0,l,... ,N-1, 



where cjj = 2, j = 0, 1, . . . ,N — 2; aN-i = 1- In this case 



ei (t, (j)j 



At' 



e2{tAj)=' 



2VN 
(0 
t^ 



A8VN 



if j = l,3,... ,N-1, 
-a,-;U,(l + (3 + /i,)2)+0(i3) if j = 0,2,... ,N-2, 

if j = l,3,... ,iV-l, 
-ajfj,j{l + 2{3 + fijf 

+ {lj] + 5fi, + 7f)Y^\o{t^) if j = 0,2,... ,N-2, 



2 7r(j + l) 
2N 



(3.8) 



(3.9) 



for j = 0,1, . . . ,N — 1, are the eigenvalues of the operator A^. 



where Hj = —4 sin 

The relevant estimates for the Euler ( ^.3| ) and KN (3^) methods are of the form 



Si{t,g) = f,f- + 0{t') 



^2 



S2{t,g) = f.',- + 0{t'). 



12 



(3.10) 
(3.11) 



Comparing the above estimates we see that, given the same order of approximation, 
the error of the OE method is much greater for eigenfunctions with small numbers and 
much smaller for higher harmonics. 

Observe that the error in the classical schemes comes from the difference between the 
eigenvalues of the step operator and exp{tAK), while their eigenfunctions coincide. The 
step operators of OE algorithm have error in both the eigenvalues and the eigenfunctionq^. 
However, as shown numerically, the eigenvalues of the operators Si{t) and 5*2 (t) approxi- 
mate the spectrum of exp{tAK) better than the eigenvalues of the classical schemes. 

To make sure the above is true, let us find spectrum Xj{t), j = 0,1, . . . , N—1, of the step 
operator 5*2 (t) for Laplacian Ad. We show in Appendix that for j odd the eigenvalues 
and eigenfunctions of operators 5*2 (t) and exp{tAu) coincide. This is exactly why for 
these harmonics the errors of the OE methods ( |3.8D and ( |3.9D vanish. The remaining -g- 



where Cji'^) fulfills the "dispersion" 



eigenvalues are X2j{t) = exp(t^j(t)), j = 0, 1, 
equation ( [A.9| ). 

The results of numerical solution of this equation are shown in Fig. 1 which provides 
the values of \^j — iJ,2j\ as function of j (curve 1). 

Clearly, for the majority of harmonics the eigenvalue error of the OE method is much 
smaller than the error of the KN scheme (curve 2). Besides, one should note the "unifor- 
mity" of the spectrum estimate of scheme ( |2.9| ): the error weakly depends on the number 
of the eigenvalues. A similar situation holds also for the operator iAu. 
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Figure 1. Error in even eigenvalues of step operator for tlie OE metliod (1) and the Krank- 
Nickolson nictlrod (2) for one-dimensional Laplace operator with Dirichlet boundary conditions 
t 



\,N = 1024. 



rel. error 
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Figure 2. Relative error of solution of the Schrodinger equation for the OE method (1) and the 

Krank-Nickolson method (2) with random initial vector, t — i, N — 128. 



This property of the step operator in the numerical scheme is important when solution of 
the input evolution problem includes contributions from all eigenfunctions of operator Ak- 



This is the case, for example, in solving problem (2.4) with skew-Hermitian operator Ak 
( S chr odinger equation ) . 



^It is interesting to observe that the spectra of Si{i) and 52 (t) coincide and their eigenfunctions differ 
only at boundary points. 
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Fig. 2 shows relative errors of the numerical solution of the problem 

^ = iAnf, f{0)=g 

as a function of time for the OE scheme ( |2.9|) (curve 1) and the Krank-Nickolson one 
(curve 2). The initial vector g is chosen as a random one; it is uniformly distributed on 
the unit sphere in C . 

4 Conclusion 

Splitting methods, including the operator exponential one, are widely used for solving 
difference linear and quasilinear initial-boundary value problems ||2|, ^. The proposed 
modification of the OE method can be applied when the evolution operator is represented 
as a sum of a difference operator with constant coefficients on a rectangular domain in Z* 
(the main part) and some, perhaps nonlinear, operator (perturbation). If the boundary 
conditions for the main part do not allow explicit computation of the input operator 
exponential, the problem can be approached by a splitting method in two stages: first we 
split off the perturbation, and then calculate an approximate exponential of the main part 
using the method proposed in this work. 

Consider a rectangular domain in Z'^ for the Laplacian with boundary conditions of 
the 3rd kind. Even in this case application of methods like the implicit Euler or Krank- 
Nickolson schemes requires iteration procedures to obtain the resolvent. The method we 
propose is explicit and, as follows from the one-dimensional examples provided in the 
work, competitive with conventional methods. 

One of important features of our method is a "uniform" property of the spectral es- 
timate of the initial problem. We have succeeded in applying the scheme described to 
solve one- and two-dimensional Ginzburg-Landau equation p^ , four-order diffusion equa- 



tion 1 11] and some other problems. 



Appendix 

In what follows we derive the equation for the eigenvalues of operator 5*2 (t) defined by rela- 
tion ( |2.9D for one-dimensional difference Laplace operator A/j with the Dirichlet boundary 
conditions. To simplify the calculations, we assume that A^ is even: A^ = 2M. 

In this case the "boundary" operator Grl is of the form: Gkl = — 2(5o, where Qq 
is the orthogonal projection on vector xq = ^^ ^~^ , where {ej} -^ is a standard basis 

in C . It is easy to calculate the exponential of such an operator: 

exp {\tGKi\ =E+ (e-* - l) Qo- (A.l) 



2 
Let A be an eigenvalue of S2{t) and ^x the corresponding eigenfunction. Then 

52(t)$A = exp i\tGKL\ exp {tko) exp (^^G^l j ^a = A^a- (A.2) 
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Using expression (AA) and notation "^x = exp i^tGKLj ^x, we express the relation in 
the form: 



A^ - exp(tAi3) Wa = A (1 - e^*) Qq^ 



X- 



(A.3) 



Let {/j}j=o ti^ ^^^ eigenbasis of operator Al (Laplacian with periodic boundary condi- 
tions): 



1 f2TTJk\ 

exp ( I — ^^^ 1 



the corresponding eigenvalues I'j being 



j,k = 0,1, 



.N-1, 



-4sm -, 



J = 0, 1, , 



,A^-1. 



(A.4) 



(A.5) 



Note that Uj = vn-j and Vj = fi2j+i for j = 1,2, .. . ,M, where //fc are the eigenvalues 
of Af). Note also that the right hand side of ( |A.3| ) for any value of ^a is proportional to 
vector xo whose expansion with respect to basis ( |A.4| ) is: 



N-l 



N-1 



Xo = ^{xo,fj)fj 



j=Q 



'2N 



i=o 



E i+»p -9^ /. 



,27r_j 

N 



(A.6) 



If Cj, where j = 0, 1, . . . , A'^ — 1, are the expansion coefficients of function ^\ with respect 



to basis (A.4), then relation (|A.3| ) can be expressed in the form: 

(A - exp(tz^j))cj = axX (l - e^*) (xq, fj), j = 0,l,... ,N - 1, 



(A.7) 



where ax = (^A,a^o)- 

In order to find all solution of equation ( A. 61 ) we consider two cases: 

1) a\ = 0. In this case Cj = or A = exp(tz>'j) for each j = 0, 1, . . . , A^ — 1. Since 

^\ 7^ 0, there exists index / such that q ^ 0. In this case A = exp(tf/) = exp(tz^7v-i) foi' 

/ = 1,2,... ,M and Cj = for j ^ l,N — I. The corresponding eigenfunctions are found 

from the condition a\ = Q: 



T /,N • 27rZ / 1 



0,1, 



A-1, 



/ = 1,2, 



M. 



(A. 



Note that these functions are the eigenfunctions of operator A/j corresponding to the 
eigenvalues /X2«-i) where I = 1,2, . . . , M. 

2) a\ 7^ 0. In this case A ^ exp{tvj) for each j = 0, 1, . . . , A — 1 as follows from ( A.7 ). 
Multiplying both sides of ( [A.6| ) by x_J^^Uy.\ and summing over j we obtain: 



N-l 



ax = axX (1 - e^*) ^ 



A — exp(tz^j 



j=0 ^^ 3' 

With condition oa 7^ and also (A.5) and (A.6) this relation takes the form: 

A/-1 



1 



1 



,2t 



M 



1 1 sp _ 

1 - exp(te) ^ 2 ^^ 1 



4 + z^,- 



exp(i(i/j - 0) / ' 



(A.9) 
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where ^ = In j. It is easy to see that equation ( A^) has exactly M real roots, one in each 
interval {uj, fj+i), where j = 0, 1, . . . , M — 1. These roots can be easily found numerically 
by the bisection method. 

Denote the solutions of ( |A.9| ) by ij{t), where j = 1, 2, . . . , M. Then, for the spectrum 
Aj(t), where j = 0, 1, . . . , A^ — 1, of operator S2{t) we finally obtain: 

\2j{t) = exp(t^j(t)), A2J+1 = exp(tz^j), i = 0, 1, . . . , M - 1. 

Note that for the eigenvalues \{t) = exp(ft^(t)) of the step operator S2{t) corresponding 
to (Schrodinger) operator iA/j equation (|A.9|) is of the form: 



tant / 1 , 1^^^ 4 + z/,- 



1 1^^^ 



M 1 tan {\ti) 2 ^ tan {\t {vj - ^)) 
The solutions of this equation are in good agreement with those of (A. 9). 
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